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^^ , This paper presents new results about the optimization based generation of chemical reac- 

tion networks (CRNs) of higher deficiency. Firstly, it is shown tliat the graph structure of the 
realization containing the maximal number of reactions is unique if the set of possible com- 

J^ ' plexes is fixed. Secondly, a mixed integer programming based numerical procedure is given for 

S . computing a realization containing the minimal/maximal number of complexes. Moreover, the 

linear inequalities corresponding to full reversibility of the CRN realization are also described. 
The theoretical results are illustrated on meaningful examples. 
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1 Introduction: chemical reaction networks and their significance 

Positive (nonnegative) systems are characterized by the property that all state variables remain 
positive (nonnegative) if the trajectories start in the positive (nonnegative) orthant. Thus, positive 
systems play an important role in fields such as chemistry, economy, population dynamics or even 
in transportation modeling where the state variables of the models are often physically constrained 
to be nonnegative [9]. It is remarked that many non-positive systems can be transformed into the 
positive class either through invertible coordinates transformations or by using other approaches, 
where the distortion of the phase-space can be kept minimal in the region of interest [26] . 

Chemical Reaction Networks (CRNs) form a wide class of positive (or nonnegative) systems 
attracting significant attention not only among chemists but in numerous other fields such as physics, 
or even pure and applied mathematics where nonlinear dynamical systems are considered [30]. 
Beside pure chemical reactions, CRNs are often used to model the dynamics of enzymatic systems 
[5], intracellular processes, metabolic or cell signalling pathways [15]. The increasing interest towards 
reaction networks among mathematicians and engineers is clearly shown by recent tutorial and 
survey papers [27, 1, 6]. 

It is known from the so-called "fundamental dogma of chemical kinetics" that reaction networks 
with different graph structures and even with different sets of complexes might generate identical 
dynamical system models (i.e. sets of differential equations). This means that CRNs with struc- 
turally very different reaction mechanisms can show exactly the same behaviour in the state space 
that is usually the space of chemical specie concentrations. However, many strong analysis results 
of chemical reaction network theory (CRNT) depend on the graph structure of the studied CRN. 
There is a clear need therefore to define and search for distinguished structures among the possi- 
ble alternatives. The integration of logical expressions into mixed integer programming problems 
[24, 4] has opened the possibility to formulate the computation of certain reaction structures with 
advantageous properties as an optimization problem [28]. 

Mixed Integer Nonlinear Programs (MINLPs) are the most general constrained optimization 
problems with a single objective. These problems can contain continuous and integer decision vari- 
ables without any limitations to the form and complexity of the objective function or the constraints. 
As it is expected, the solution of these problems is rather challenging [12]. A special subset of op- 
timization problems is the class of Mixed Integer Linear Programs (MILPs) where the objective 
function and the constraints are linear functions of the decision variables. Effective solvers have 
been developed for MILPs, although it is known that their solution is NP-hard. In the chemical and 



biochemical fields, efficient combinatorial optimization algorithms are widely applied e.g. in perma- 
nental polynomial computation [19], metabolic pathway construction, control analysis or metabolic 
network reconstruction [3]. It is noted that the evolutionary approach can also be very successful 
in solving complex chemically originated optimization problems [18]. 

In [28], the notion of realization was introduced for the unique definition of a reaction net- 
work, and a mixed integer linear programming (MILP)-based numerical procedure was proposed to 
compute sparse and dense realizations of mass-action reaction networks corresponding to the same 
mathematical model, solving important part of a problem that was raised almost 30 years ago in 
[16]. The purpose of this paper is to present new results in the field of optimization based generation 
of reaction network structures. 

2 Basic notions and tools 

2.1 Structural and dynamic description of CRNs obeying the mass action law 

The overview in this subsection is largely based on [28]. A CRN obeying the mass action law is a 
closed system under isothermal and isobaric conditions, where chemical species Xj, i = 1, ...,n take 
part in r chemical reactions. The concentrations of the species denoted by Xi, {i = l,...,n) form 
the state vector, i.e. Xi = [Xj]. The elementary reaction steps have the following form: 

n n 

^Oij-Xj -> ^ftjXj, j = l,...,r (1) 

4 = 1 j=l 

where aij is the so-called stoichiometric coefficient of component X, in the jth reaction, and /3j£ is 
the stoichiometric coefficient of the product X^. The linear combinations of the species in eq. (1), 
namely Yl^=i (^ij^i ^^d Y17=i Pij^i ^^ J — 1; • • • j'^ ^^'^ called the complexes and are denoted by 
Ci, C2, . . . , Cm- Note that the stoichiometric coefficients are always nonnegative integers in classical 
reaction kinetic systems. The reaction rates of the individual reactions can be described as 

n n 

P,=k,Yl[X,r^=k,llx^^^ , J = l,...,r (2) 

i=i 1=1 

where kj > is the reaction rate constant of the jth reaction. 

If the reactions Ci — t- Cj and Cj — >■ Ci take place at the same time in a reaction network for 
some i,j then this pair of reactions is called a reversible reaction (but it will be treated as two 
separate elementary reactions). 



Similarly to [10], we can assign the following directed graph (see, e.g. [2]) to the reaction network 
(1) in a straightforward way. The directed graph D = {Vd,Ed) of a reaction network consists of 
a finite nonempty set V^ of vertices and a finite set E^ of ordered pairs of distinct vertices called 
directed edges. The vertices correspond to the complexes, i.e. Vd = {Ci, C2, . . . C^}, while the 
directed edges represent the reactions, i.e. {Ci,Cj) G Ed if complex Cj is transformed to Cj in the 
reaction network. The reaction rate coefficients kj for j = 1, ... ,r in (2) are assigned as positive 
weights to the corresponding directed edges in the graph. Where it is more convenient, the notation 
fc^- will be used for denoting the reaction rate coefficient corresponding to the reaction Cj -^ Cj. 

A set of complexes {Ci, C2, ■ ■ ■ ,Ck} is a linkage class of a reaction network if the complexes of 
the set are linked to each other in the reaction graph but not to any other complex [11]. There are 
several possibilities to represent the dynamic equations of mass action systems (see, e.g. [10], [14], 
or [8]). The most advantageous form for our purposes is the one that is used e.g. in Lecture 4 of 
[10], i.e. 

± = Y-Ak- ^(x) (3) 

where x G M" is the concentration vector of the species, Y € M"^™- stores the stoichiometric 
composition of the complexes, A^ € M"^^™ contains the information corresponding to the weighted 
directed graph of the reaction network, and ^ : M" 1— > M™ is a monomial-type vector mapping 

defined by 

n 

ijj{x) = llxf\ j = l,...,m (4) 

i=l 

where yij = \Y]ij. It is remarked that the numerical solution of kinetic differential equations can 
be a challenging task requiring advanced integration approaches [29]. The exact structure of Y and 
Ak is the following. The ith column of Y contains the composition of complex Cj, i.e. Yji is the 
stoichiometric coefficient of Cj corresponding to the specie Xj. Ak is a column conservation matrix 
(i.e. the sum of the elements in each column is zero) defined as 

In other words, the diagonal elements [Akju contain the negative sum of the weights of the edges 
starting from the node Cj, while the off-diagonal elements [j4fc]jj, i ^ j contain the weights of the 
directed edges {Cj,Ci) coming into Cj. Based on the above properties, it is appropriate to call A^ 
the Kirchhoff matrix of a reaction network. 



To handle the exchange of materials between the environment and the reaction network, the so- 
called "zero-complex" can be introduced and used which is a special complex where all stoichiometric 
coefficients are zero i.e., it is represented by a zero vector in the Y matrix (for the details, see, e.g. 
[10] or [7]). 

We can associate an n-dimensional vector with each reaction in the following way. For the 
reaction Cj — )> Cj, the corresponding reaction vector denoted by /i^ is given by 

h, = [y].,, - [YU (6) 

where \Y].^i denotes the ith column of Y. Similarly to reaction rate coefficients, whenever it is more 
practical, h'-- denotes the reaction vector corresponding to the reaction Cj -^ Cj. 

The rank of a reaction network denoted by s is defined as the rank of the vector set H = 
{hi,h2 ■ ■ ■ ,hr} where r is the number of reactions. The elements of H span the so-called sto- 
ichiometric subspace denoted by S, i.e. S = span{/ii, /12 . . . , /ir}. The positive stoichiometric 
compatibility class containing a concentration xq is the following set [11]: 

[xQ + s)nRl 

where M^ denotes the positive orthant in M^. 

The deficiency d of a reaction network is defined as [10, 11] 

d = m — I — s (7) 

where m is the number of complexes in the network, / is the number of linkage classes and s is the 
rank of the reaction network. 

A reaction network is called reversible, if each of its reactions is a reversible reaction. A reaction 
network is called weakly reversible, if each complex in the reaction graph lies on at least one directed 
cycle (i.e. if complex Cj is reachable from complex Cj on a directed path in the reaction graph, 
then Cj is reachable from Cj on a directed path). An important point of the well-known Deficiency 
Zero Theorem [11] says that the ODEs of a weakly reversible deficiency zero CRN are globally 
stable with a known logarithmic Lyapunov function for all positive values of the reaction rate 
coefficients. Therefore (among other realization problems) it is of interest whether we can find a 
(weakly) reversible deficiency zero kinetic realization of a nonnegative polynomial system. 

Using the notation M = Y ■ A^, eq. (3) can be written in the compact form 

± = M ■ il}{x) (8) 

The invariance of the nonnegative orthant for CRN dynamics is shown e.g. in [6]. 



2.2 Kinetic realizability of positive (nonnegative) polynomial systems 

An autonomous polynomial nonlinear system of the form 

X = fix) (9) 

is called kinetically realizable or simply kinetic, if a mass action reaction mechanism given by eq. (3) 
can be associated to it that exactly realizes its dynamics, i.e. f{x) =Y-Ak- ip{x) where ip contains 
the monomials, matrix Y has nonnegative integer elements and A^ is a valid Kirchhoff matrix (see 
section 2.1 for its properties). In such a case, the pair (Y,Ak) will be called a realization of the 
system (8) (note that Y contains all information about the composition of the monomials in -0 in 
the case of mass-action dynamics). As it is expectable from linear algebra, the same polynomial 
system may have many parametrically and/or structurally different realizations. Thus, two CRNs 
will be called dynamically equivalent if they realize the same polynomial system of the form (9). 
Therefore, CRN A will also be called a realization of CRN B,\i A and B are dynamically equivalent. 
The problem of kinetic realizability of polynomial vector fields was first examined and solved in 
[16] where the constructive proof contains a realization algorithm that produces the directed graph 
of a possible associated mass action mechanism. It is important to remark here that the above 
mentioned realization algorithm typically produces high deficiency CRNs that are non-minimal in 
the sense that they usually contain more reactions and complexes than the minimal numbers that 
are necessary to realize the given kinetic polynomial system. According to [16], the necessary and 
sufficient condition for kinetic realizability is that all coordinates functions /j of the right hand side 
of (9) must have the form 

fi{x) = -Xigi{x) + hi{x), i = l,...,n (10) 

where gi and hi are polynomials with nonnegative coefficients. 

2.3 Mixed integer linear programming 

A special subset of optimization problems is the class of Mixed Integer Linear Programs (MILPs) 
where the objective function and the constraints are linear functions of the decision variables. A 
mixed integer linear program with k variables (denoted hy w ^M.^) and p constraints can be written 



as [22]: 



... T 

inmimize c w 



subject to: 

Aiw = hi 

A2W < 62 (11) 

h ^Wi < Ui for i = 1 , . . . , A; 

Wj is integer for j (z I, I ^ {1, . . . ,k} 

where c G M^, Ai G MPi><*^, A2 G M^^xfe^ ^^^^^ ^^ +^3 = P- 

If all the variables can be real, then (11) is a simple linear programming problem that can be 
solved in polynomial time. However, if any of the variables is integer, then the problem becomes 
NP-hard. In spite of this, there exist a number of free (e.g. YALMIP or the GNU Linear Program- 
ming Kit) and commercial (such as CPLEX or TOMLAB) solvers that can efficiently handle many 
practical problems [17, 20, 21]. 

A propositional logic problem, where a statement denoted by S must be proved to be true given 
a set of compound statements containing so-called literals Si, ... , 5"^, can be solved by means of a 
linear integer program. For this, logical variables denoted by 5i {5i G {0, 1}) must be associated 
with the literals Si. Then the original compound statements can be translated to linear inequalities 
involving the logical variables 6i [25, 4]. 

2.4 Computing CRN realizations with the minimal/maximal number of reac- 
tions as a MILP problem 

For convenience, this subsection briefly summarizes the results of [28] without going into the details. 
The starting point is that a kinetic polynomial system of the form (8) is given with its parameters. 
This means that M is known, the stoichiometrix matrix Y is also known from the monomials of ip, 
and we would like to determine the Kirchhoff matrix A^ G M™^"^ that fulfils given requirements. 
The characteristics of the mass-action dynamics can be expressed in the form of the following 



equality and inequality constraints: 

Y-Ak = M (12) 

m 

J^[^fc],j=0, j = l,...,m (13) 

[Ak\ij > 0, i, j = 1, . . . , m, i / j (14) 

[^fc]n<0, i=l,...,m (15) 

where the decision variables are the elements of Aj.. Clearly, constraints (13)-(15) express that 
we are searching for a valid Kirchhoff connection matrix. To make the forthcoming optimization 
problems computationally tractable, appropriate upper and lower bounds are introduced for the 
elements of A]f_ as 

< [Ak]ij <lij, i, i = 1, . . . , m, i / j (16) 

ki<[Ak]ii <0, i = l,...,m. (17) 

In this problem set, we are searching for such Af^ that contains the minimal /maximal number of 
nonzero off-diagonal elements. For this, we introduce logical variables denoted by 5 and construct 
the following compound statements 

6ij = 1 o [Aklij > e, i,j = l,...,m, i y^ j (18) 

where the symbol "o" represents "if and only if", and < e <C 1 (i.e. elements of A), below e are 
treated as zero). Taking into consideration (16), statement (18) can be translated to the following 
linear inequalities (see, e.g. [4]) 

< [Aklij - eSij, i,j = 1,... ,m, i j^ j (19) 

< -[Ak]ij + kjSij, i,j = 1,... ,m, i^ j (20) 

Now we are able to compute the realization containing the minimal/maximal number of reactions 
by minimizing/maximizing the objective function 

m 

Ci(5)= Y. ^^0 (21) 

«,i = 1 

The realizations of a reaction network containing the minimal and maximal number of reactions 
will be called the sparse and dense realizations, respectively [28]. 



2.5 A simple motivating example 

Consider the simple reaction mechanism depicted in Fig. la). It is easy to check that the reaction 
structures in Figs 1 b) , c) , d) and e) lead to the same dynamical description as the original structure 
a), namely 

xi = 3kiX2 — k2Xi 

X2 = -3kixl + k2xl, (22) 

with ki,k2 > 0, 5/c2 > ^i (i-e. the CRNs in Fig. 1 are dynamically equivalent). It is worth having a 
look at the structural properties of the different realizations of eq. (22) shown in the subfigures. The 
realizations in Figs, l.a) and b) are irreversible, the structure in Fig. l.c) is weakly reversible, while 
the networks in Figs, l.d) and e) are fully reversible. The deficiencies of the first four realizations 
a)-d) are 1, while the deficiency of realization e) is zero. This means that both the weaker Deficiency 
one theorem and the stronger Deficiency zero theorem can be applied to all realizations a)-e), and 
this way to the dynamical system described by eq. (22) (see [11]). Shortly speaking, the Deficiency 
one theorem for such weakly reversible networks as c) says that its differential equations admit 
precisely one steady state in each positive stoichiometric compatibility class. Moreover, by applying 
the Deficiency zero theorem to realization e), we obtain the additional valuable fact that each steady 
state of (22) is asymptotically stable within the corresponding positive stoichiometric compatibility 
class with the Lyapunov function: 

y(x) = ^x.(ln(|i)-l)+x:, (23) 

where x* denotes the equilibrium point of (22) corresponding to the given stoichiometric compati- 
bility class. 

First of all, the above example shows very transparently that important structural properties 
such as deficiency, reversibility or weak reversibility are not encoded uniquely in the polynomial 
differential equations of a kinetic system. Secondly, it is definitely of interest to develop computa- 
tional tools to search for realizations with such properties that are useful in the dynamical analysis 
of given kinetic polynomial systems or CRNs. 
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Figure 1: Dynamically equivalent reaction networks 

3 Properties of dense realizations and its consequences 

The main result of this section is that the dense realization of a CRN is structurally unique if the 
set of possible complexes is fixed. We recall that the realizations of a reaction network containing 
the minimal and maximal number of nonzero reaction rate coefficients are called the sparse and 
dense realizations. 



3.1 The uniqueness of the structure of dense realizations 

Firstly, we state the following result. 

Theorem 3.1. // a set of kinetic differential equations denoted by S is given with matrices M and 
Y , then the directed unweighted graph of any realization of S must he a subgraph of the directed 
unweighted graph of the dense realization. 

Proof. The proof is based on the following elementary fact of linear algebra. Consider an inhomoge- 
nous set of linear equations: 



Ax = b 



(24) 



10 



If p is any specific solution of (24) then the entire solution set of (24) can be given as 

{p + V \ V is any solution of Ax = 0} (25) 

The matrix equation Y ■ A^ = M (see eqs. (3) and (8)) obviously defines m sets of linear equations 
of the form 

Y■[Ak].,^ = [M].,„ i = l,...,m (26) 

where the unknown is [^fc]-,i (i-e-, the ith column of A^)- For any i, let us assume that p = [Ak].^i 
is a dense solution of (26) i.e., it contains the maximal possible number of nonzero elements. Let 
I denote the number of nonzero elements in p. If / = tti then the theorem is trivial, so from now 
on we assume that I < m. Let us assume furthermore that p' ^ p is also a solution of (26) and let 
(ji, . . . ,jq) denote the indices where p{jk) = while p'{jk) 7^ for /c = 1, ... g. With Q' > 0, this 
means that the directed unweighted reaction graph defined by p' is not a subgraph of the directed 
unweighted reaction graph defined by p. Then, according to (25), p' can be written as 

p' = p + v, (27) 

where Y ■ v = {). According to our assumption, w 7^ 0, and necessarily, v{ji) ^ 0, . . . ,v{jq) 7^ 0. 
Let /' denote the number of nonzeros in p'. Since /' < / (because p is a dense solution), there must 
exist indices {hi, . . . ,hz) disjoint from (ji, . . . ,jq) with z > q such that v{hk) = —p{h].) 7^ for 
k = 1, . . . ,z. Then for any X & M., p" = p + X ■ v is also a solution of (26), and A can always be 
chosen such that p" contains more nonzero elements than p, which is clearly a contradiction. D 

We note that we did not use the further restriction that [^fc]-,i is an appropriate column of a 
Kirchhoff matrix, but this was not needed for the proof. Now we easily obtain our following result 
about the uniqueness of the dense realization. 

Theorem 3.2. // a set of kinetic differential equations denoted by T, is given with matrices M and 
Y , then the directed graph structure of its dense realization is unique. 

Proof. The proof is the special case of the proof of Theorem 3.1 with /' = / and q = z. D 

3.2 Important consequences and special cases 

The following remarks contain some important immediate consequences and additions to Theorems 
3.1 and 3.2. 

11 



Rl According to Theorems 3.1 and 3.2, dense realizations give a unique "superstructure" for a 
CRN in the sense that the reactions of any reahzation of a CRN must form a subset of the 
reactions of the dense reahzation if the set of possible complexes is given. In other words, 
reactions that are not present in the dense realization cannot appear in any other realization. 

R2 Obviously, dense realizations are parametrically not unique. There may exist several dense 
realizations for a CRN with different reaction rate constants (weights) but always with the 
same graph structure. 

R3 The graph structure of a CRN with a given set of complexes is unique if and only if the graph 
structures of its sparse and dense realizations are identical. 

This fact is easy to see: If the structures of the dense and sparse realizations are identical, 
then it directly follows that the graph structure of the CRNs is unique, since the only possible 
unique structure is determined by the dense realization (that is the sparse realization at the 
same time). In other words, any realization of the CRN can contain neither more nor less 
reactions than the dense realization does, the structure of which is unique. If the graph 
structure of the CRN is unique, then it trivially implies that the structures of the dense and 
sparse realizations are identical. 

R4 The dense realization of a CRN is not only a theoretical construction but it can be practically 
determined using well-formulated numerical procedures that are treatable even in the case of 
several hundred complexes and species (see, e.g. [13, 25]). 

R5 Sparse realizations of CRNs are structurally not unique, there may exist several sparse real- 
izations for a given CRN with different graph structures (see later in subsection 5.1). 

4 Transforming additional constraints corresponding to preferred 
CRN properties into linear inequalities 

This section presents some further answers to the open problems originally set in [16] from an 
optimization point of view. 
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4.1 Computing realizations with the minimal/maximal number of complexes 

In this section, the detailed MILP formahsm will be presented for computing CRN realizations that 
contain the minimal/maximal number of complexes from a predefined complex set. 

Let us assume again that the set of feasible complexes is a'priori given with matrix Y. The 
constraints written in eqs. (12)-(17) corresponding to the characteristics of mass-action dynamics 
are used here again without change. Then, the minimization or maximization of the number of non- 
isolated complexes in the reaction graph is based on the following simple observation. A complex 
disappears from the reaction network's graph, if both the corresponding column and row in A^ 
contain only zeros. This means that no directed edges start from or point to this complex in the 
graph and therefore it becomes an isolated vertex that can be omitted. 

For the optimization, m boolean variables denoted hy 5i, i = 1, ... ,m are introduced. Using 
these boolean variables, the following compound statements are introduced: 

m m 

5i = l^ Y^ [Ak\ij,+ Y. [M2,i>0, i = l...,m (28) 

jl = l J2 = l 

Eq. (28) means that the value of bi is 1 if and only if there is at least incoming/outgoing directed 
edge in the reaction graph to/from the ith complex. For practical computations, the statement (28) 
is modified as follows: 

ra m 

5^ = l^ Y [Mi,h + Y [^fc]j2,i > e, i = 1 . . . , m (29) 

ii=i J2=i 

where again < e <C 1 (see, eq. (18)). Using the bound constraints (16)-(17), the linear inequalities 
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corresponding to (29) are the following 

m rn 



(30) 






J2 = l 
J2 #i 



/ 



0<- J2 i^kkn- J2 [^fc]j2,» + e- 



/ . ''ijl + 2^ ^2 4 ~ 

ii=i J2=i 

y ji#'' j2#i y 

Now, the objective function to be minimized or maximized can be written as 






J2 = l 



\ 



Si, i = 1, . . . , m (31) 



C^2(5) = X^5. 



(32) 



4 = 1 



In contrast to the algorithm summarized in section 2.4, minimizing/maximizing the number of 
non- isolated complexes is not straightforward to parallelize (see also [28]). However, the number 
of integer variables in this case is only m, compared to m'^ — m when minimizing/maximizing the 
number of reactions. 

4.2 Computing reversible realizations 

Here, the basic constraints (12)-(17) expressing the properties of mass action dynamics and lower 
and upper bounds for the reaction rate coefficients will be used again for the optimization. To 
distinguish between zero and nonzero reaction rate coefficients, a small positive scalar e is applied 
again, similarly to the previous case described in section 4.1. 

The additional constraint for the full reversibility of the CRN structure is not difficult to for- 
mulate as 

[Ak]ij > £2 o [Ak]j,i > €2, Vi > J. (33) 

where e2 is a positive threshold value such that e < e2- The linear inequalities equivalent to (33) 
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can be written as 



< (62 - e) - [Ak]^, + {kj - €2) • 4'^ , Vi > j (34) 

< (62 - e) - [Ak]ji + (Iji - €2) • 5f^ , Vi > j (35) 

Q<[Ak\ij-e2-5f^, yi>j (36) 

0<[Aky-e2-dff, yi>j (37) 



where lij is the upper bound for [j4fc]jj as it is introduced in eq. (16). Furthermore, IHi!^ — l integer 
variables are introduced for the representation of the reversibihty constraint that are denoted by 

In order to exclude reaction rate coefficients between e and e2, and to obtain a numerically stable 
solution, the following additional constraints in the form of a compound statement are introduced 

[Ak]ij < e OR [Ak]ij > £2 + 7, (38) 

where 7 is a small positive threshold value that is in the same order of magnitude as £2- The set of 
inequalities equivalent to (38) is given by 

0<4f^ i^J (39) 

< kj - [Akhj - ihj - e) ■ 6§\ i^j (40) 

0<[Ak]i,-{e2+j)-6ff, i^j (41) 

0<-6^+sl^+6i^, i^j (42) 

0<sif-6lf, i + 3 (43) 

0<5lf-<5j;\ i + 3 (44) 

where j'^^, b^^' and b^^' represent altogether 3(m^ — ?n) integer variables. 

It is remarked that the inequalities (34) - (37) and (39) - (44) express only constraints and no 
objective function is associated to reversibility in itself. However, the reversibility constraints can 
be easily combined with the minimization/maximization of either the number of reactions or that of 
the non-isolated complexes, still in the framework of mixed integer linear programming. Moreover, 
the strict reversibility constraint can be modified into the minimization/maximization of reversible 
reactions in a straightforward way. It is emphasized finally, that the constraints presented in this 
subsection together with an appropriate MILP solver are suitable for deciding whether a reversible 
realization exists for a given CRN or not. 

15 



5 Examples 

For the examples described in this section, the YALMIP modehng tool was applied under the MAT- 
LAB computational environment [20] using both the freely available GLPK [21] and the commercial 
CPLEX [17] solvers. 

5.1 Non-uniqueness of sparse realizations 

For the illustration of several structurally different sparse realizations of a reaction network, let us 
recall a literature example that was originally published in [8]. The original CRN with all reaction 
rate coefficients equal to 1 is shown in Fig. 2 a). The CRNs in Figs. 2 b) and c) were obtained by 
using the parallel and non-parallel version of the method described in [28] and summarized in section 
2.4, respectively, using the GLPK MILP solver. The network shown in Fig. 2 d) was computed using 
the CPLEX solver, using a non-parallel approach. These results show that additional constraints 
in the optimization procedure may be used to select the required sparse realization from the set of 
possible alternatives. 



5.2 Motivating example continued 

Consider again the reaction network shown in Fig. 1 a) with parameters ki 
matrices characterizing the CRN realization are the following. 
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M = Y -A, 



3 
-3 



-2 
2 



(46) 



5.2.1 Computing a realization with the minimal number of complexes 

Finding a realization with the minimal number of complexes using the method described in section 
4.1 with parameters lij = 100 Vi,j and e = 10~® gives the following result: 
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Figure 2: Original reaction network and its three different sparse realizations. Only reaction rates 
different from 1 are indicated. 

(2) (2) 

It's straightforward to check that M = Y ■ Aj^ . Here, Aj^ gives a deficiency structure that is 
shown in Fig. 1 e). 
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5.2.2 Computing a dense reversible realization 

If we search for a reversible realization given by eq. 46 that contains the maximal number of 
nonzero reaction rate coefficients (i.e. a dense reversible realization), we have to combine constraints 
(12)-(17), (34)-(37), (39)-(44) and (19)-(20), and maximize the objective function (21). Using the 
parameters e = 10"'^, €2 = 0.05, 7 = 0.01 we obtain a fully reversible structure given by the following 
Kirchhoff matrix 



A 



(3) 



-1.0200 0.6467 33.3333 
0.9600 -0.7067 66.6667 
0.0600 0.0600 -100.0000 



(48) 



which gives a deficiency 1 structure shown in Fig. l.d. Again, it's clear that Y ■ A^ = Y ■ A 



(3) 



5.3 Equivalent reversible realization of an irreversible reaction network 

Let us start from the reaction network that is depicted in Fig. 3. This network contains 9 complexes, 
2 linkage classes and 8 irreversible reaction steps. The rank of the stoichiometric subspace is 3, 
therefore the deficiency of the network is 4. The matrices characterizing the network are given by 
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(49) 



(50) 



Running the algorithm described in section 4.2 with parameters e 



10-^ 



£2 = 0.05, 7 = 0.01, 



where the objective function to be minimized was the number of nonzero reaction rate coefficients. 
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Figure 3: Irreversible reaction network with a deficiency of 4 



gave tlie following Kirchhoff matrix: 
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It is again easy to verify that 



Y-Ak = Y-A'k 







0.5 
-3.5 

1.5 



0.5 

0.5 

-0.5 









(51) 



(52) 



The above result implies that the deficiency zero theorem can be applied to the dynamics of the 
original irreversible reaction network shown in Fig. 3. Moreover, due to the existence of a defi- 
ciency reversible realization with linearly independent reaction-pairs, the dynamics of the reaction 
networks exhibit a dissipative Hamiltonian structure as it was shown in [23]. 

6 Conclusions 

Different possible realizations of dynamically equivalent CRNs have been studied in this paper 
with the help of mixed integer linear programming. The main contributions of the paper can be 
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Figure 4: Zero deficiency reversible reaction network dynamically equivalent to the one shown in 
Fig. 3 

summarized as follows. Firstly, it has been shown that the structure of a so-called dense realization 
of a given CRN is unique, and the structure of any other realization is the subgraph of the dense 
realization if the set of complexes is given. By computing a possible sparse realization, it is also 
possible to test numerically, whether the structure of a CRN is unique or not. Secondly, a method 
has been given for finding a CRN realization with the minimal number of complexes (from within 
a predefined set) in the framework of MILP. Finally, the numerically feasible constraints (linear 
(in)equalities) for determining reversible realizations of CRNs have been presented. The theoretical 
findings have been illustrated on examples. The results clearly show the power of linear programming 
combined with propositional logic for determining preferred realizations of reaction kinetic systems. 
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